Coho Sampling Map with Leaflet

Building an interactive map of Coho sampling sites and Conservation Units using the leaflet R package

Load libraries:

library(leaflet)
library(tidyverse)
library(leaflet.extras)
library(htmltools)
library(htmlwidgets)
library(readr)

Read in the file with site names, CU designation, and coordinates:

sites <- read_tsv("./CohoSamplingMapCU/data/ALL_SITES_UPDATEDwCOLOURS.txt", col_names = TRUE)
head(sites)
## # A tibble: 6 × 9
##   SITE  FULL       REGION   CU    LATITUDE LONGITUDE pop_cu     color       n
##   <chr> <chr>      <chr>    <chr>    <dbl>     <dbl> <chr>      <chr>   <dbl>
## 1 AAL   Aaltanhash BC       CO-30     53.1     -128. AAL(CO-30) #FDCDAC    18
## 2 ALB   Albreda    Thompson CO-09     52.5     -119. ALB(CO-9)  #A65628    34
## 3 ALO   Alouette   BC       CO-47     49.3     -123. ALO(CO-47) #FB8072    45
## 4 ASH   Ashlu      BC       CO-10     49.9     -123. ASH(CO-10) #B3CDE3    16
## 5 AVO   Avola      Thompson CO-09     51.8     -119. AVO(CO-9)  #A65628    42
## 6 BAR   Barriere   Thompson CO-09     51.2     -120. BAR(CO-9)  #A65628    38

Plot Points

Let's start with a map of sampling locations with two interactive features:

  1. a label with site code when you hover over the points
  2. a popup with site code, region, and sample size when you click on a point
sites %>% 
  leaflet() %>%
  addProviderTiles("CartoDB") %>% 
  addCircleMarkers(popup = ~paste0("<b>", SITE, "</b>", "<br/>", "Region: ", REGION, "; n = ", n), label = ~SITE, radius = 5, color = 'black', fillOpacity = 1, stroke = T, weight = 0.8, fillColor = 'tomato')
## Assuming "LONGITUDE" and "LATITUDE" are longitude and latitude, respectively

addProviderTiles() loads a base map from a number of different providers available in the leaflet package. Use names(providers) to see a list.

names(providers)[1:20]
##  [1] "OpenStreetMap"               "OpenStreetMap.Mapnik"       
##  [3] "OpenStreetMap.DE"            "OpenStreetMap.CH"           
##  [5] "OpenStreetMap.France"        "OpenStreetMap.HOT"          
##  [7] "OpenStreetMap.BZH"           "OpenSeaMap"                 
##  [9] "OpenPtMap"                   "OpenTopoMap"                
## [11] "OpenRailwayMap"              "OpenFireMap"                
## [13] "SafeCast"                    "Thunderforest"              
## [15] "Thunderforest.OpenCycleMap"  "Thunderforest.Transport"    
## [17] "Thunderforest.TransportDark" "Thunderforest.SpinalMap"    
## [19] "Thunderforest.Landscape"     "Thunderforest.Outdoors"

Points are plotted using addCircleMarkers() and the interactive features associated with the points are included here as arguments, with ~ to point to columns in sites :

  • popup: popup text that appears on click
  • label: text box that appears when you hover the cursor over a point

Other arguments are included to modify the aesthetics of the points.

Next, I want to colour the points according to the conservation unit that each site is currently managed in, and add 'CU' to the popup. I included a colour column in the sites dataframe; this is a list of colours for each CU that I used for consistency across all figures in the paper. We'll use that column to create a color palette associated with CU and apply to the map.

pal <- colorFactor(palette = unique(sites$color), levels = unique(sites$CU))

sites %>% 
  leaflet() %>%
  addProviderTiles("CartoDB") %>% 
  addCircleMarkers(popup = ~paste0("<b>", SITE, "</b>", "<br/>", "Region: ", REGION, "<br/>", "CU: ", CU, "<br/>", "n: ", n, "<br/>"), label = ~SITE, radius = 6, fillColor = ~pal(CU), fillOpacity = 0.8, stroke = T, color = "black", weight = 0.8)
## Assuming "LONGITUDE" and "LATITUDE" are longitude and latitude, respectively

Add Polygons

Now I'll add the CU polygons, which are stored in a .shp file. A couple of spatial packages are needed to read and manipulate this file.

library(rgeos)
library(rgdal)

cu.map <- readOGR("./CohoSamplingMapCU/data/Coho_Salmon_CU_simplified/cu.simp.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Volumes/Storage/epic4/EPIC4/MAPS/CohoSamplingMapCU/data/Coho_Salmon_CU_simplified/cu.simp.shp", layer: "cu.simp"
## with 44 features
## It has 1 fields
## Rows: 44 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): Full_CU_IN, SPECIES_QU, CU_name, Type
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

This is a SpatialPolygonsDataFrame and the metadata is stored in the data slot.

head(cu.map@data)
## # A tibble: 6 × 4
##   Full_CU_IN SPECIES_QU CU_name        Type         
##   <chr>      <chr>      <chr>          <chr>        
## 1 CO-01      CO         Boundary Bay   VREQ[Current]
## 2 CO-04      CO         Lillooet       Current      
## 3 CO-05      CO         Fraser Canyon  Current      
## 4 CO-07      CO         Lower Thompson Current      
## 5 CO-08      CO         South Thompson Current      
## 6 CO-09      CO         North Thompson Current

I created another colour palette for the CUs that's separate from the points. There are some CUs that were unsampled; I included these boundaries on the map as well, shaded white and with a black border.

color.list <- read.table("./CohoSamplingMapCU/data/color_list.txt", header = F, stringsAsFactors = F)  # list of unique colours for CUs
head(color.list)
##        V1
## 1 #B2DF8A
## 2 #CAB2D6
## 3 #E5D8BD
## 4 #D95F02
## 5 #6A3D9A
## 6 #A65628
# create a second list with "#FFFFFF" (white) switched to "black"; one list will be for the polygon fill, the other for the border.
color.listb <- color.list
color.listb$V1[which(color.listb$V1 == "#FFFFFF")] <- 'black'

# fill palette associated with CU names stored in the data slot of the SpatialPolygonsDataFrame
cu.pal <- colorFactor(palette = color.list$V1, levels = cu.map@data$Full_CU_IN)

# border palette associated with CU names stored in the data slot of the SpatialPolygonsDataFrame
cu.palb <- colorFactor(palette = color.listb$V1, levels = cu.map@data$Full_CU_IN)

Let's add the CU polygons to the map with sampling points using addPolygons(). We need to specify the dataset (in this case, the SpatialPolygonsDataFrame called cu.map). The highlightOptions lets us send the polygons to the back so we can still see the labels and popups for the points.

sites %>% 
  leaflet() %>%
  addProviderTiles("CartoDB") %>% 
  addCircleMarkers(popup = ~paste0("<b>", SITE, "</b>", "<br/>", "Region: ", REGION, "<br/>", "CU: ", CU, "<br/>", "n: ", n, "<br/>"), label = ~SITE, radius = 6, fillColor = ~pal(CU), fillOpacity = 1, stroke = T, color = "black", weight = 0.8) %>%

addPolygons(data = cu.map, weight = 0.5, 
              fillColor = ~cu.pal(Full_CU_IN), 
              fillOpacity = 0.3, 
              color = ~cu.palb(Full_CU_IN),
              highlight = highlightOptions(sendToBack = T)
              )
## Assuming "LONGITUDE" and "LATITUDE" are longitude and latitude, respectively

Let's add a label and popup for the CUs and have some fun with the highlightOptions. I want to have the polygons stand out when I move the mouse over them.

sites %>%
  leaflet() %>%
  addProviderTiles("CartoDB") %>%
  addCircleMarkers( popup = ~paste0("<b>", SITE, "</b>", "<br/>", "Region: ", REGION, "<br/>", "CU: ", CU, "<br/>", "n: ", n, "<br/>"), label = ~SITE, radius = 6, fillColor = ~pal(CU), fillOpacity = 1, stroke = T, color = "black", weight = 0.8) %>%
  
  addPolygons(data = cu.map, weight = 0.5, 
              fillColor = ~cu.pal(Full_CU_IN), 
              fillOpacity = 0.3, 
              color = ~cu.palb(Full_CU_IN),
              popup = ~Full_CU_IN,
              highlight = highlightOptions(weight = 1, 
                                           color = "black",
                                           sendToBack = T,
                                           bringToFront = F, opacity = 1,
                                           fillOpacity = 0.2)
              ) 
## Assuming "LONGITUDE" and "LATITUDE" are longitude and latitude, respectively

Extra Map Features

There are tons interactive features in the leaflet package to make the map even more fun. For more information: https://rstudio.github.io/leaflet/

I'm just going to add a few more things:

  • a default map extent
  • a button to reset the map to the default extent (use it to go back after zooming in)
  • a search box (magnifying glass): use the site code to identify it on the map and zoom to that location (for example, try searching SAL, YAK, SNO). Typing one letter into the search box a drop down list should appear)
  • a layers control check box to toggle off/on the CU boundaries and to show only the sites in either of the two main regional groups (BC and Thompson); this will be at the top right corner. Note that I added group = ~REGION inside addCircleMarkers() and group = ~CU inside addPolygons() to define the layer groups.
sites %>%
  leaflet() %>%
  addProviderTiles("CartoDB") %>%
  fitBounds(lng1 = -135, lat1 = 47, lng2 = -118, lat2 = 58) %>%
  addResetMapButton() %>%
  addCircleMarkers( popup = ~paste0("<b>", SITE, "</b>", "<br/>", "Region: ", REGION, "<br/>", "CU: ", CU, "<br/>", "n: ", n, "<br/>"), label = ~SITE, radius = 6, fillColor = ~pal(CU), fillOpacity = 1, stroke = T, color = "black", weight = 0.8, group = ~REGION) %>%
  
  addSearchFeatures(targetGroups = ~REGION, options = searchFeaturesOptions(zoom = 8, openPopup = TRUE)) %>%
  
  addPolygons(data = cu.map, weight = 0.5, 
              fillColor = ~cu.pal(Full_CU_IN), 
              fillOpacity = 0.3, 
              color = ~cu.palb(Full_CU_IN),
              popup = ~Full_CU_IN,
              highlight = highlightOptions(weight = 1, 
                                           color = "black",
                                           sendToBack = T,
                                           bringToFront = F, opacity = 1,
                                           fillOpacity = 0.2),
              group = "CU") %>%
  addLayersControl(overlayGroups = c("BC", "Thompson", "CU"))
## Assuming "LONGITUDE" and "LATITUDE" are longitude and latitude, respectively